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BLOCK TRIANGULAR PRECONDITIONING FOR STOCHASTIC 

GALERKIN METHOD 

BIN ZHENG, GUANG LIN, AND JINCHAO XU 

Abstract. In this paper we study fast iterative solvers for the large sparse linear systems 
resulting from the stochastic Galerkin discretization of stochastic partial differential equations. A 
block triangular preconditioner is introduced and applied to the Krylov subspace methods, including 
the generalized minimum residual method and the generalized preconditioned conjugate gradient 
method. This preconditioner utilizes the special structures of the stochastic Galerkin matrices to 
achieve high efficiency. Spectral bounds for the preconditioned matrix are provided for convergence 
analysis. The preconditioner system can be solved approximately by geometric multigrid V-cycle. 
Numerical results indicate that the block triangular preconditioner has better performance than the 
traditional block diagonal preconditioner for stochastic problems with large variance. 
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l/"~) 1. Introduction. Stochastic partial differential equations (SPDEs) provide math- 

ematical models for uncertainty quantification in many physical and engineering ap- 

^* plications, including flows in heterogeneous porous media [18], thermo-fluid processes 

ySL [2S1 [U|, flow-structure interactions [3D] , etc. These models are proposed to quantify 

the propagation of uncertainties in the input data, such as coefficients, forcing terms, 

t£^ boundary conditions, initial conditions, and geometry of the domain, to the response 

quantities of interests. 

Numerical methods for solving SPDEs can be categorized as either non-intrusive 
or intrusive methods. Non-intrusive methods are based on sampling techniques where 
a number of uncoupled deterministic partial differential equations need to be solved. 
Examples of such methods are Monte Carlo method and stochastic collocation method. 
Monte Carlo method is the most straightforward approach and its convergence is in- 

iy\ dependent of the number of stochastic dimensions. It is useful for problems with very 

£■ — , large stochastic dimensions where most of the other methods suffer from the "curse of 

dimensionality" . The disadvantage of Monte Carlo method is that it has a very slow 

^+ rate of convergence, i.e., proportional to 1/yN where N is the number of samples. 

Stochastic collocation method combines a Galerkin approximation in the physical 
space and a collocation/interpolation in the zeros of suitable orthogonal polynomi- 
als in the probability space. It achieves fast convergence for problems with moder- 
t> ate stochastic dimensions and smooth solutions in the stochastic domain 33.P221I5]. 

k> However, the number of collocation points in a tensor grid grows exponentially with 

respect to the number of random variables. To overcome this problem, the Smolyak 
sparse grid stochastic collocation method is developed for problems with high random 
dimensions |38j . The main advantage of the aforementioned non-intrusive methods 
is the ease of implementation, namely, existing deterministic solvers can be employed 
without any modification to solve the uncoupled deterministic problems. 

Stochastic Galerkin method |15) is considered as an intrusive method in the sense 
that it results in coupled systems which cannot be solved by deterministic solvers 
directly. One advantage of the stochastic Galerkin method is that the number of 
equations is relatively small and scales as w 1/2 P (p is the order of the stochastic dis- 
cretization) times the number of sparse grid stochastic collocation equations [38l [41] , 
which is appealing for high order stochastic discretization. In [111 [5]. experimental 
comparisons of the stochastic Galerkin method and stochastic collocation method are 
provided which show that the stochastic Galerkin method is advantageous in terms of 
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computational cost when efficient solvers are available. Stochastic Galerkin method 
is based on the concept of generalized polynomial chaos expansion J37J 131] which 
provides an exponentially convergent approximation for modeling smooth stochastic 
processes. It applies Galerkin projection onto a space of generalized polynomial chaos 
in the probability space. Such stochastic discretization is often combined with a finite 
element discretization in the physical space. The resulting method is often called 
the stochastic Galerkin finite element method for SPDEs. However, the size of the 
coupled linear system grows dramatically when increasing the solution resolution in 
both stochastic space and physical space. Hence, the development of fast solvers is 
necessary for solving these coupled linear systems efficiently. Fortunately, the matrix 
of this coupled linear system has particular structures which can be utilized in de- 
signing fast solvers. In particular, although deterministic solvers can not be applied 
directly to the coupled linear system, they may be used as building blocks for the 
design of new solvers. 

There are a number of studies on fast solvers for the large coupled linear system 
including preconditioned Krylov subspace methods and multigrid methods. We re- 
fer to [551 GH] for comprehensive overviews and comparisons of iterative solvers for 
stochastic Galerkin discretizations. Due to the sparsity pattern of the linear system, 
Krylov subspace methods that require only matrix- vector multiplication are very at- 
tractive. Preconditioning techniques have been studied to accelerate the convergence 
of the Krylov subspace methods. The block diagonal preconditioner (also known as 
mean-based preconditioner) for the conjugate gradient (CG) method is the most con- 
venient one that has been observed to be robust for problems with small variance 
[Til [55"] . In [55] , theoretical eigenvalue bounds for the block diagonal preconditioned 
system matrix are derived. Also, in [35], a symmetric positive definite Kronecker 
product preconditioner has been introduced which makes use of the entire informa- 
tion contained in the coupled linear system to achieve better performance in terms 
of the CG iteration counts. Multigrid methods with optimal order of computational 
complexity in the physical space are also investigated theoretically and numerically 
for stochastic Galerkin discretizations [02 GUI EH [57] . The hierarchical structure in 
the stochastic dimension resulting from the use of hierarchical polynomial chaos basis 
functions is explored in the design of iterative solvers [21 [53J [55J • 

In this work, we combine the optimality of the multigrid method in physical space 
with the hierarchical structure in probability space to obtain efficient preconditioner 
for the Krylov type iterative methods. In particular, we propose a block triangu- 
lar preconditioner for the generalized minimal residual (GMRes) method, and the 
generalized preconditioned conjugate gradient (GPCG) method. We also consider a 
symmetric block preconditioner for the standard conjugate gradient (CG) method. 
Numerical results indicate that block triangular preconditioner is more efficient than 
the traditional block-diagonal preconditioner for problems with large random fluctu- 
ations. 

The rest of the paper is organized as follows. We start in Sect ion [2] from a descrip- 
tion of a model elliptic problem with random diffusion coefficient and an overview of 
the corresponding stochastic Galerkin finite element discretization. In Section [3] we 
summarize the structures of the stochastic Galerkin matrix. In Section HI we intro- 
duce the block triangular preconditioner and its symmetrized version. Convergence 
analysis of the preconditioned Krylov subspace methods is described in Section [5] 
Finally, the performance of the proposed block triangular preconditioner for GMRes 
and GPCG method are demonstrated in Section [U 



2. Model problem and stochastic Galerkin discretization. We consider 
the following elliptic problem 

- V ■ (a(x, u)Vu(x, w)) = f(x), x e D c R d , (2.1) 

u(x,u>) = 0, x G dD, 

where the diffusion coefficient a is a real- valued random field defined on D, i.e. for each 
x G D, a(x, •) is a random variable with respect to the probability space (17, IF, P). 
We assume that a is bounded and uniformly coercive, i.e. 



^min i "max 



G (0,+oo) : P(u G Q : a(x,u) G [a min ,a max ], Vx e D) = 1. (2.2) 



Introducing the tensor product Hilbert space V = L 2 P (Q.) ® Hq(D) with inner 
product defined by 



(u, v)y= / I / Vu(x,oj) • Vw(x,cj) da; I dP(w), 
the weak solution u G V is a random function such that V v G V: 



a(x, w)V«(x, u>) ■ Vv(x,uj) dx I dP(w) = / I / f(x)v(x, u>)dx I dP(w). 

(2.3) 



The well-posedness of the above variational problem (2.3) follows from (2.2) and the 
Lax-Milgram lemma. 

2.1. Karhunen-Loeve expansion. The first step of the stochastic discretiza- 
tion is to approximate the input random field a(x, ui) by the truncated Karhunen- 
Loeve expansion 



a(x,to) sa a m (x,oj) := a(x) + ^ \^kb k (x)^ k (uj), (2.4) 



fc=i 



where a(x) is the mean value of a(x,uj), X k and b k (x) are the eigenvalues and eigen- 
functions of the integral operator C : L 2 (D) — > L 2 (D) defined by J D Cov a (a;, •) u(x) dx. 
Given mean value a(x) and a continuous covariancc function Cov a (a;, y), the truncated 



Karhunen-Loeve expansion (2.4) approximates a(x,uj) with minimized mean square 
error |15j . The number of terms in the expansion, m, is determined by the eigenvalue 
decay rate, and in turn, depends on the stochastic regularity, i.e., the smoothness of 
the covariance function. The expansion coefficients Cfe( w ) are pairwise uncorrelated 
random variables with images T k = 6c (^): and probability density functions (PDFs) 
Pk : Tfe —> W l . The joint PDF of the random vector £ = (£i, . ..,£m) is denoted by 
p(0, and the image F = H-f =l T k . 

If a(x,u>) is a Gaussian random field, £& will be independent Gaussian random 
variables, and the joint PDF />(£) = H™ =1 Pk{£,k)- hi general, for non-Gaussian random 
field, £/c are not necessarily independent and their distributions are not known. Several 
methods have been developed to estimate the distributions of £fc and to simulate non- 
Gaussian processes using Karhunen-Loeve expansion, see |24[ 136] . A non-Gaussian 
random field may also be approximated by polynomial chaos expansion, see |13| 126]. 



An example of the covariance function is given by the exponential covariance function 

Cov a (x,y)=a 2 exp(~\x~y\/L), (2.5) 
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where a is the standard deviation and L is the correlation length. 

Remark 1. When replacing the diffusion coefficient a(x,u>) by the truncated 
Karhunen-Loeve expansion a m (x,uj), it is important to verify the uniform coercivity 



condition (2.2) so that the problem 



— V • (a m (x, lo)Vu(x, w)) = f(x), x 6 D, (2-6) 

u(x,oj) — 0, x G 3D. 

is well-posed. For more discussions, including the estimate of the error between the 



two solutions of (2.1) and (2.6), we refer to \ 121 [21/. 

One advantage of using Karhunen-Loeve expansion is the separation of the stochas- 
tic and deterministic variables for the stochastic function a{x,uj). In addition, by the 



Doob-Dynkin lemma [17] . the solution of (2.6) can be written in terms of £, i.e. 
u(x,u>) = u(x,£,i(lu), . . . ,^ m (uj)). The stochastic problem (2.6) is then reformulated 
as the following deterministic parametrized problem: 

- V • (a m (x, OVufo 0) = /(*), x€D,£<=T, (2.7) 

u{x,0 = 0, xe 3D, £e r. 

2.2. Stochastic Galerkin discretization. Since the weak solution u(x,£) is 
defined in a tensor product space V, we consider finite dimensional approximation 
space also in tensor product form, i.e. Vh, p = Xh <£> £ p . When the solution is 
smooth/analytic in stochastic variables, spectral approximation using global poly- 
nomials of total degree < p in m variables defined in T 

E p = span{^i(0, ■ • • , <M0} C L 2 p (T) 

are good candidates for approximations in the stochastic space. The global polynomi- 
als ipi are chosen to be orthogonal polynomials associated with the density function 
p, often referred to as generalized Polynomial Chaos (gPC) [3S], e.g., Legendre poly- 
nomials for uniform distribution, Hermite polynomial for Gaussian distribution, etc. 
The dimension of the space S p is given by the following formula 

{m + p)\ 

iv ? = TT - - 

m\p\ 

For example, when m = 6, p = 4, dim(S4) = 210. 

For the spatial approximation, we choose the standard finite element space, i.e., 
space of piecewise polynomials with respect to a given mesh Dh (h is the spatial 
discretization parameter) 

X h = span{0i(a:),...,0;v !t (a;)} C Hq(D), 

where N x is the dimension of Xh- Hence, the discrete solution Uh, P can be written as 
the following polynomial chaos expansion 



«/.,p(z.o = 51)^(^(0 - E Z)^>.(*) ^(0- (2.8) 

j=l j=l Vs=l / 

An a priori error estimate for the error \\u— Uh, P \\ l 2 (r)^H 1 (d) is given in [3]. Statistical 
information including mean, variance, etc., can then be obtained from the explicit 
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formula given by Eqn. (2.8), which give good approximations to those of the exact 
solution. 

The stochastic Galerkin finite element method is obtained by applying the Galerkin 
projection in the tensor product space Vh <p - More precisely, find Uh, P G Vh, p such that 



where 



and 



B(uh, P ,v) 



B(uh, P , v) = (/, v), V v e Vh- 



P(0 / a m (x,C)^xU h , P (x,0 ■ \7 x v(x,£)dxd£, 



(f,v) 



P(0 / v(x,$f(x)dxdZ. 

r Jd 



Since B(-,-) is a symmetric and positive definite bilinear form, it also introduces an 
inner product and the associated norm denoted by (•, -) a and || • || a , respectively. 
We refer to [3] for a more thorough discussion of the stochastic Galerkin method. 

3. Matrix structures and iterative solvers. In order to design efficient and 
robust iterative solvers, it is important to study the structure of the corresponding 
matrix. For example, the Galerkin matrix A is symmetric and positive definite (SPD) 
by the uniform elliptic assumption (2.2). As a consequence, the CG method can be 
applied to solve the linear system. Since the solution space V and the approximation 
space Vh are both tensor spaces, the matrix A also contains a tensor product structure. 
Moreover, A has block sparsity structure and hierarchical structure as described below. 

3.1. Tensor product structure. Consider the semi-discretization in stochastic 
domain, the corresponding Galerkin projection u^(-,£) : T — » X^ satisfies for each 



stochastic basis polynomial ipi(£), i = 1, 



,N f 



-V-(o TO (a ! ,OV«W(a!,0)^(Op(Ode= ffr)MtMt)*$, (3-1 



-V ; 



where the semi-discrete approximation u( p '(x,t;) = Ylj=i u j( x )' l Pj(0- Note Eqn. (3.1 1 



,Ne 



is a system of N$ equations with N^ unknown functions {uj(x)}^ ly i.e. the stochastic 
'stiffness' matrix is given by 



A = 



( M t i 



Ax.2 
A 2 ,2 



Ax 
A 



N e 



\ 



2.N e 



\ An ( ,i AN i ,2 ■■■ An^N,: J 
where each entry Ai.j contains spatial differentiation and is given by 



A, 



-V ■ (a m (x,^u 3 (x))^ 3 (0M0p(0^. 



Substitute in the truncated Karhunen-Loeve expansion (2.4) for a m , we get 
Ai tj = -V • (a(x)Vuj(x)) J ifc(0lfc(0p(0 d £ 



III n 

fc=l Jr 



Next, we apply Galerkin projection to discretize spatial differential operators. Let 



(x) =J2 U J,sMx), 



u , - 

8=1 



where {<f) s (x)} s ^ 1 is the set of finite element basis functions. Multiply the spatial 
derivative terms in each Ai.j by the test function 4> r (x) (for r = 1, . . . ,N X ) and 
integration by parts gives 

y^U^ s \ a{x)V4> s (x) ■ V(j> r (x)dx, S~] y/^k^U^s / b k (x)V(j) s (x) -V<t) r (x)dx. 

> JD i- < s =l J° 

(3.2) 
From (3.2 1, we define spatial stiffness matrices Kq and Kk (for k = 1, . . . , m) as 

i^o( r > s) = / a(x)V()> s (x) ■ V(f> r (x) dx, 
Jd 



k=l 



Kk{r, s) = VA/c / b k (x)\7<f) s (x) ■ \/<p r {x) dx, r,s = l,..., N x . (3.3) 

JD 

Similarly, we define the stochastic matrices Go and {Gk} k n =1 as 

G (i,j) = f 1>M)MtMS)4t, 

G k (i,j) = j ^j(0^i(Op(0^, i,j = l,...,N ( . (3.4) 

Finally, the Galerkin matrix A can be written in terms of Kronecker products: 

i n 
A = G ®K Q + ^G k ®K k . (3.5) 

fc=i 

In practice, one does not assemble the Galerkin matrix A explicitly. Instead, using 
the tensor product structure (3.5), only to + 1 spatial stiffness matrices Kk and m + 1 
stochastic matrices Gk need to be stored. 

3.2. Block sparsity structure. By construction, the gPC basis functions {4>i}i Jj 
are orthonormal with respect to the PDF p, i.e. 

'r 

Furthermore, the stochastic matrices Gk (k = 0, 1, . .. ,m) are sparse as can be seen 
from the explicit formulas for the matrix elements (3.4) given in [28] (Theorem 2.2). 
As a consequence, the Galerkin matrix A has a particular block sparsity structure, 
see Fig |3.1| for the case when p = 4, m = 4. This block sparsity structure is essential 
for the construction of many efficient iterative solvers for stochastic Galerkin methods 
as seen in [231 122] ■ 

3.3. Hierarchical structure. In addition to the block sparsity property, the 
matrix A also has a hierarchical structure due to the hierarchical nature of the stochas- 
tic approximation spaces S p , i.e., S p _i C S p , and the use of gPC basis functions which 
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Fig. 3.1. Block sparsity structure of the Galerkin matrix A (m ■■ 



4j. 



are also hierarchical. Namely, the basis functions of higher order stochastic approx- 
imation space are obtained by adding homogeneous high order polynomials while 
keeping basis functions from lower order space. We remark that the hierarchical 
basis functions are also used in many finite element computations for deterministic 
problems, see [5] and the references therein. 

Let A be the Galerkin matrix corresponding to X/j <£> S p with gPC basis of order 
p. It can be rewritten as 



A = 



A W T 
W D 



, A = Kq, 



(3.6) 



where A is the Galerkin matrix corresponding to .Xfc®S p _i, W represents the coupling 
of higher-order and lower-order stochastic modes, and D is a block diagonal matrix 
corresponding to the homogeneous gPC basis of degree p. 

In Q31 [S3] , several hierarchical 



This hierarchical structure is shown in Fig 3.1 



approaches have been discussed assuming weak coupling of different orders of ap- 
proximation. Multilevel methods in the stochastic domain based on this hierarchical 
structure have been designed and used as preconditioners for Krylov methods in [2 8) 
which exhibited good convergence properties. 

4. Block preconditioners. Preconditioning techniques are necessary in order 
to accelerate the convergence of the Krylov subspace methods when applied to ill- 
conditioned linear systems. Since the Galerkin matrix A has a block sparse structure, 
it is natural to consider block preconditioners. 

4.1. Block-diagonal preconditioner. It is known that the block-diagonal pre- 
conditioner (also known as the mean-based preconditioner) defined by 



B, 



Go ® K , 



(4.1) 



works very well with the CG method when the variance of the diffusion coefficient 
a(x, uS) is small P31[53]. Spectral bounds for the block-diagonal preconditioned system 
matrix have been derived in 125 . 



By the uniform ellipticity assumption, the mean stiffness matrix Kq is SPD. 
Multigrid method can be used to invert each diagonal block approximately. 

4.2. Block-triangular preconditioner. Another choice is to use the lower 
block triangular part of the matrix A p as a preconditioner. This can be viewed as 
applying one step of the block Gauss-Seidel method with zero initial guess. Consider 
a splitting of the stochastic approximation space S p , 

and the corresponding splitting of the global approximation space 
V h , p = (X h <g> S p _i) 8 (X h <g> (3p\S p _i)), 



which results in a 2 x 2 block structure given by ( |3.6| . 

Given u^ ' = 0, for k = 1, 2, . . . , the block Gauss-Seidel iterate u^ k+1 ^ is given by 
the following two steps: 

• Find u( fc +V2) G u (k) + x h (g> S p _i such that 



)= / P(£) / v(x,£)f(x)dxde, Vv&X h ®Z p _ 



B{u( k+l ' 2 \ 
• Find u( fc+1 ) e u( fc +V 2 ) + E h ® (S p \3 p _ 1 ) such that 

B(u< fc >, «)= /p(0 / «(»,0/(»)dxd$, Wel^fH^.!). 

In matrix notation, the above block Gauss-Seidel method can be described by the 
following matrix splitting 



A = 



A 
W D 



-W T 




We define the block triangular preconditioner 
The corresponding preconditioner system 



A 
W D 



(4.2) 



A 
WD 




. u 2 _ 


= 


Fi 

. F 2 . 



may be solved inexactly by the standard multigrid V-cycle. 

Remark 2. The block triangular preconditioner Bj- may also be motivated by 
considering the block L U factorization 



A = 



A 







' I A- l W T 


W 


S D 




/ 



, S D = D-WA~ 1 W 1 



(4.3) 



It is known that with the "ideal" block triangular preconditioner 



Bt — 



A 

w S D 



the GMRes method converges in at most two iterations. However Bt is impractical 
because the Schur complement Sjy is computationally expensive to invert. Replacing 
Sd by D in Bt results in the block triangular preconditioner Bt ■ 
Remark 3. A different form of the block LU factorization 



A = 



I W T D- 1 
I 



Sa 






D 



I 
D- l W I 



S A = A- W T D- 1 W 1 (4.4) 



is studied in [32] from which a hierarchical Schur complement preconditioner is de- 
rived. 

Since Bt is nonsymmetric, we can use it with the GMRes method [25] or GPCG 
method [7j. To apply the standard PCG method, we may consider the block symmet- 
ric Gauss-Seidel method as the preconditioner, i.e. 



B 8 := 



' A " 




' A " 


-i 


WD 




D 





A W T 
D 



It is clear that B$ is SPD and the standard PCG method is guaranteed to converge. 

5. Convergence Analysis. In this section, we give eigenvalue bounds for the 
matrix preconditioned by block triangular preconditioner which is crucial to the con- 
vergence of the preconditioned Krylov subspace methods. We note that the spectral 
properties of block triangular preconditioner for saddle point problems have been 
studied in QUEHO]. 

The following two lemmas are useful. 

Lemma 5.1. \31^ The eigenvalues of AB^f are positive real numbers, and the 
spectrum satisfies 

a{AB T x ) C {1}U o-{S,D) 

where S — D — W A~ l W T is the Schur complement of A in A, and o~(S, D) contains 
the eigenvalues \x corresponding to the generalized eigenvalue problem 



Sz = \jlDz. 



(5.1) 



Lemma 5.2. [25] Let K Q and Kf- be the stiffness matrices defined in \3.3ty , then 
the eigenvalues of Kq K^ belong to the interval 



fc«fe 



a*^ 



-y^kWbkWoc, -V^fcll^fcl 
a a 



(5.2) 



depending on the positivity of bk (x) . 

Next, we give the main result of this section. 

Theorem 5.3. The spectrum of the preconditioned Galerkin matrix AB^f satis- 
fies 

o(AB T l ) C (0,1]. 



Proof. By Lemma 5.1 it suffices to show that o~{S, D) C (0, 1] 
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Let \i be an eigenvalue of the generalized eigenvalue problem (5.1), and z be the 
corresponding eigenvector, 



(D - WA- 1 W T )z = f iDz 

=* z T (D- WA~ 1 W T )z = pz T Dz 

z T (D-WA- l W T )z 
z 1 Dz 

Since both S and D are symmetric positive definite, we know /j, > 0. 



(5.3) 



Similarly, from (5.3), 



(1- fj,)Dz= WA- 1 W T z 
=*> (l-fj,)z T Dz = {W T z) T A- 1 (W T z) 
(W T z) T A-\W T z) 



=>■ 1 — jX = 



z T Dz 



and the symmetric positive definiteness of D and A , we conclude that 

1 - /i > =>> /i < 1. 

D 

For the case p = 1, the following result gives a better lower bound for fi. 

Theorem 5.4. When p — 1, i/ie eigenvalue /i of the preconditioned Galerkin 
matrix AB^, satisfies 



>:L - c Ei2 A * 



k lloo ' 



fe=l 

Proo/. Notice ^(ABy 1 ) = e^P^A), and 

-l 



B^A 



A 
VK D 



A W T 
W D 



A 
W D 



W T 




Let /x £ a(B T 1 A), v — [fijVa] be the corresponding eigenvector. We have 

(m-i) 



W^ T 




Vl 

V-2 



A 
W D 



t'2 



which implies 



Hence, 



Notice that 



(A- W l D- x W)v x = nAvi. 



(I - A-'W 1 D~ L W)vx = fivi. 



A = K , W = Y^G k ®K kl D = G ®K 



k=l 
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where Ko, Kk are defined by (3.3 1, Gk is the first column of Gk defined by (3.4), and 
Go 



Imxm assuming the stochastic basis functions are normalized. 
Using the properties of the Kronecker product, we get 

m m 

D~ l W = Y, G k ® K^Kk, A- X W T = Y, Gj ® i^ 1 ^, 



fe=i 



/=i 



and 



i^W^L^PF = ^] 



1=1 

m m 



(Gf ® ^o" 1 ^) IJ2 Gk ® K o lRk 
\fc=i 

n 

i=i fe=i 

m 

^(G^G fe ) ® (i^o- 1 ^) 2 



fe=i 



where c is a constant depending on the stochastic basis functions. 
Applying Lemma [572] to conclude that 



H>l-c 



m 1 



2 

oo ' 



fe=l 



D 

For p > 1, we do not have a similar estimate for the lower bound of ijl due to the 
lack of simple formulae for A^ 1 which is the inverse of a sum of Kronecker products. 
Furthermore, it is well known that the eigenvalue information alone is not sufficient 
to predict the convergence behavior of the GMRes method [8]. On the other hand, 
we observe the uniform convergence rate with respect to both mesh size h and the 
stochastic discretization parameters p and m in the numerical experiments described 
in Section |6] 

6. Numerical results. In this section, we evaluate the performance of the block 
triangular preconditioner and compare with the block diagonal preconditioner. The 
test problem is taken from [5S]. Let D = (-0.5,0.5) x (-0.5,0.5). The covari- 
ance function is given by (2.5) where a = 0.1, 0.2, 0.3, or 0.4 and L = 1. The 



independent random variables £j are assumed to be uniformly distributed with range 
7t = ( — V3) v3)j for i — I, ■ ■ ■ ,m. The mean value of the random coefficient a(x) = 1, 
and the right hand side f(x, y) = 2(0.5 — x 2 — y 2 ). 

In our implementation of the stochastic Galerkin method, the linear finite element 
method is used for spatial discretization on the uniform triangulation with mesh size 
h. Multivariate Legendre polynomials of total degree < p are used for the stochastic 
discretization. For each Krylov subspace method, the outer iteration starts from a 
zero initial guess and ends when the relative residual error in Eucliean norm less than 
10~ 10 is achieved. All computations are done in MATLAB on a laptop with a 2.8 
GHz Intel processor and 4 GB of memory. 

In the following tables, we use the abbreviations, B-GS : block Gauss-Seidel 
method; BD-PCG: PCG with block diagonal preconditioner; BT-GPCG: GPCG[1] 
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Block-GS 
^-BD-PCG 
-*-BT-GPCG 
-B-BT-GMRes 
-l-BS-PCG 




Fig. 6.1. Number of iterations with exact solve of diagonal blocks (p = 4, m = 6, h = 1/64J. 

method (for details see [7]) with block triangular preconditioner; BT-GMRes: GM- 
Res(lO) method (restarted every 10 iterations if needed) with block triangular pre- 
conditioner; BS-PCG: PCG with block symmetric Gauss-Seidel preconditioner. 

The sparsity information for the relevant matrices used in the block diagonal 
preconditioner Bd and the block triangular preconditioner Bt is reported in Table 

Table 6.1 
Size and number of nonzeros of the relevant matrices (p = 4, m = 6, h = 1/64J. 





A 


A 


W 


D 


B D 


size 


833,490 


333, 396 


(500, 094) x (333, 396) 


500, 094 


833, 490 


nnz 


23,863,938 


8,228,948 


6,583,136 


2,468,718 


4,114,530 



In Table [672] we report the performance of the block preconditioners with respect 
to the spatial mesh size h, and the stochastic discretization parameters p and m 
when a = 0.1 is fixed. From Table 6.2 it is shown that with the block diagonal, 
block triangular preconditioners, or block symmetric Gauss-Seidel preconditioner, the 
Krylov subspace methods converge uniformly with respect to h, p, and m. 



Table 6.2 
Number of iterations with diagonal blocks solved by one multigrid V(2, 2) cycle (a 



0.1, rn- 



or 6). 





BD-PCG 


BT-GPCG / BS-PCG 


BT-GMRes 


h 


p = 2 


P = 3 


p = 4 


p = 2 


P = 3 


p = 4 


p = 2 


P = 3 


p = 4 


1/32 


13 


13 


13 


9 


9 


9 


8 


8 


8 


1/64 


13 


13 


13 


9 


9 


9 


8 


8 


8 


1/128 


13 


13 


13 


9 


9 


9 


8 


8 


8 



Figure \6A\ shows the robustness of the different iterative methods with respect to 
the diffusivity variance a when the diagonal blocks are solved exactly. It can be seen 
that the number of iterations increases when a increases for all the iterative solvers 
considered. This seems to be related to the fact that large variance indicates strong 
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coupling between the stochastic modes of different order. Moreover, the matrix A may 
become indefinite for large a. The methods using block triangular preconditioner are 
more robust compared with the block diagonal preconditioned CG method. However, 
the cost associated with block triangular preconditioner is larger than that of the 
block diagonal preconditioner. Hence, the gain in terms of computational efficiency 



(CPU time) is not as significant as seen from iteration counts, see Table 6.3 



Table 6.3 
Number of iterations and GPU time (in seconds) with one V(2, 2) for diagonal blocks (p = 4, 
: 6, h = 1/64J. 



a 


0.1 


0.2 


0.3 


0.4 


Block-GS 


13(15.53) 


16(18.87) 


24(30.67) 


65(76.35) 


BD-PCG 


13(12.18) 


18(15.79) 


27(24.13) 


49(43.40) 


BT-GPCG 


9(12.91) 


10(13.42) 


13(17.53) 


22(30.15) 


BT-GMRes 


8(16.18) 


9(17.01) 


12(23.24) 


20(37.88) 


BS-PCG 


9(19.12) 


10(20.96) 


12(27.17) 


20(44.24) 



Table 6.3 shows the iteration counts and CPU time (in seconds) comparison of 



the five iterative solvers. Note that here all the diagonal blocks are solved inexactly by 
applying a single geometric multigrid V-cycle with two pre- and two post-smoothing 
(point Gauss-Seidel smoother). From this table we can see that block triangular 
preconditioner performs better than block diagonal preconditioner when a is large. 
We also point out that when a is small, the block Gauss-Seidel method also gives fairly 
good results comparable to those from the preconditioned Krylov subspace methods. 
The block symmetric Gauss-Seidel preconditioned CG method is computationally 
more expensive than other methods although its number of iterations is small. 

7. Conclusions. In this work we study block preconditioners for the coupled 
linear systems resulting from the stochastic Galerkin discretizations of the elliptic 
stochastic problem. The proposed block triangular preconditioner utilizes the block 
sparsity and hierarchical structure of the stochastic Galerkin matrix. The precon- 
ditioner is solved inexactly by geometric multigrid V-cycle and applied to Krylov 
subspace methods such as GMRes or GPCG. A symmetrized version of the precon- 
ditioner is proposed and applied to the standard PCG method. Numerical results 
indicate that the block triangular preconditioner achieves better efficiency compared 
to the traditional block diagonal preconditioner especially for problems with large 
variance. We also give theoretical bounds for the spectrum of the preconditioned 
system. 
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